home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / md_test.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  134 lines
  1. ;$Id: md_test.pro,v 1.6 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       MD_TEST
  8. ;
  9. ; PURPOSE:
  10. ;       This function tests the hypothesis that a sample population
  11. ;       is random against the hypothesis that it is not random. The 
  12. ;       result is a two-element vector containing the nearly-normal
  13. ;       test statistic Z and the one-tailed probability of obtaining
  14. ;       a value of Z or greater. This test is often refered to as the
  15. ;       Median Delta Test.
  16. ;
  17. ; CATEGORY:
  18. ;       Statistics.
  19. ;
  20. ; CALLING SEQUENCE:
  21. ;       Result = MD_test(X)
  22. ;
  23. ; INPUTS:
  24. ;       X:    An n-element vector of type integer, float or double.
  25. ;
  26. ; KEYWORD PARAMETERS:
  27. ;   ABOVE:    Use this keyword to specify a named variable which returns
  28. ;             the number of sample population values greater than the
  29. ;             median of X.
  30. ;
  31. ;   BELOW:    Use this keyword to specify a named variable which returns
  32. ;             the number of sample population values less than the
  33. ;             median of X.
  34. ;
  35. ;     MDC:    Use this keyword to specify a named variable which returns
  36. ;             the number of Median Delta Clusters; sequencial values of
  37. ;             X above and below the median.
  38. ;
  39. ; EXAMPLE:
  40. ;       Define a sample population.
  41. ;         x = [2.00,  0.90, -1.44, -0.88, -0.24,  0.83, -0.84, -0.74,  0.99, $
  42. ;             -0.82, -0.59, -1.88, -1.96,  0.77, -1.89, -0.56, -0.62, -0.36, $
  43. ;             -1.01, -1.36]
  44. ;
  45. ;       Test the hypothesis that X represents a random population against the
  46. ;       hypothesis that it does not represent a random population at the 0.05
  47. ;       significance level.
  48. ;         result = md_test(x, mdc = mdc)
  49. ;
  50. ;       The result should be the 2-element vector:
  51. ;         [0.459468, 0.322949]
  52. ;
  53. ;       The computed probability (0.322949) is greater than the 0.05
  54. ;       significance level and therefore we do not reject the hypothesis 
  55. ;       that X represents a random population. 
  56. ;
  57. ; PROCEDURE:
  58. ;       MD_TEST computes the nonparametric Median Delta Test. Each sample
  59. ;       in the population is tagged with either a 1 or a 0 depending on
  60. ;       whether it is above or below the population median. Samples that
  61. ;       are equal to the median are removed and the population size is 
  62. ;       corresponding reduced. This test is an extension of the Runs Test 
  63. ;       for Randomness. 
  64. ;
  65. ; REFERENCE:
  66. ;       APPLIED NONPARAMETRIC STATISTICAL METHODS
  67. ;       Peter Sprent
  68. ;       ISBN 0-412-30610-7
  69. ;
  70. ; MODIFICATION HISTORY:
  71. ;       Written by:  GGS, RSI, November 1994
  72. ;       Modified:    GGS, RSI, November 1996
  73. ;                    Fixed an error in the computation of the median with
  74. ;                    even-length input data. See EVEN keyword to MEDIAN.
  75. ;-
  76.  
  77. function md_test, x, above = above, below = below, mdc = mdc
  78.  
  79.   on_error, 2
  80.  
  81.   nx = n_elements(x)
  82.  
  83.   if nx le 10 then message, $
  84.     'Not defined for input vectors of 10 or fewer elements.'
  85.  
  86.   if nx mod 2 eq 0 then $ ;x is of even length.
  87.     ;Average Kth and K-1st medians.
  88.     medx = median(x, /EVEN) else medx = median(x)
  89.  
  90.   ;Eliminate values of x equal to median(x).
  91.   xx = x[where(x ne medx)]
  92.  
  93.   ;Values of xx above and below median(x).
  94.   ia = where(xx gt medx, above)
  95.   ib = where(xx lt medx, below)
  96.  
  97.   ;Values above are tagged 1.
  98.   if above ne 0 then xx[ia] = 1
  99.  
  100.   ;Values below are tagged 0.
  101.   if below ne 0 then xx[ib] = 0
  102.  
  103.   ;sxx = above + below
  104.   sxx = size(xx)
  105.  
  106.   h0 = where(xx eq 0, n0)
  107.   if n0 ne 0 then hi = where(h0+1 ne shift(h0, -1), nn0) $
  108.   else nn0 = 0L
  109.  
  110.   ;above median
  111.   n1 = sxx[1] - n0
  112.  
  113.   if xx[sxx[1]-1] ne xx[0] then nn1 = nn0 $
  114.   else if xx[0] eq 1 then nn1 = nn0 + 1 $
  115.   else nn1 = nn0 - 1
  116.  
  117.   if n0 eq 0 or n1 eq 0 then message, $
  118.     'x is a sequence of identically distributed data.'
  119.  
  120.   ;Formulate the Median Delta test statistic.
  121.   mdc = nn0 + nn1
  122.   if sxx[2] eq 5 then n0 = n0 + 0.0d
  123.   e = 2.0 * n0 * n1 / (n0 + n1) + 1.0
  124.   v = 2.0 * n0 * n1 * (2.0 * n1 * n0 - n0 - n1) / $
  125.       ((n0 + n1 - 1.0) * (n1 + n0)^2)
  126.   z = (mdc - e) / sqrt(v)
  127.  
  128.   ;Probability from a Gaussian Distribution.
  129.   prob = 1.0 - gauss_pdf(abs(z))
  130.  
  131.   return, [z, prob]
  132.  
  133. end
  134.